From XML data
From XML dataΒΆ
Example on how to load your data if it is currently in xml format
import json
import dolfin
from fenics_plotly import plot
import pulse
try:
from dolfin_adjoint import Constant, DirichletBC, Function, Mesh, interpolate
except ImportError:
from dolfin import Mesh, Function, Constant, DirichletBC, interpolate
# Mesh
mesh = Mesh(pulse.utils.mpi_comm_world(), "data/mesh.xml")
# Marker functions
facet_function = dolfin.MeshFunction("size_t", mesh, "data/facet_function.xml")
marker_functions = pulse.MarkerFunctions(ffun=facet_function)
# Markers
with open("data/markers.json", "r") as f:
markers = json.load(f)
# Fiber
fiber_element = dolfin.VectorElement(
family="Quadrature",
cell=mesh.ufl_cell(),
degree=4,
quad_scheme="default",
)
fiber_space = dolfin.FunctionSpace(mesh, fiber_element)
fiber = Function(fiber_space, "data/fiber.xml")
microstructure = pulse.Microstructure(f0=fiber)
# Create the geometry
geometry = pulse.HeartGeometry(
mesh=mesh,
markers=markers,
marker_functions=marker_functions,
microstructure=microstructure,
)
activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
activation.assign(Constant(0.0))
matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
activation=activation,
parameters=matparams,
f0=geometry.f0,
)
# LV Pressure
lvp = Constant(0.0)
lv_marker = markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]
# Add spring term at the base with stiffness 1.0 kPa/cm^2
base_spring = 1.0
robin_bc = [pulse.RobinBC(value=Constant(base_spring), marker=markers["BASE"][0])]
# Fix the basal plane in the longitudinal direction
# 0 in V.sub(0) refers to x-direction, which is the longitudinal direction
def fix_basal_plane(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
bc = DirichletBC(V.sub(0), Constant(0.0), geometry.ffun, markers["BASE"][0])
return bc
dirichlet_bc = [fix_basal_plane]
# You can also use a built in function for this
# from functools import partial
# dirichlet_bc = partial(pulse.mechanicsproblem.dirichlet_fix_base_directional,
# ffun=geometry.ffun,
# marker=geometry.markers["BASE"][0])
# Collect boundary conditions
bcs = pulse.BoundaryConditions(
dirichlet=dirichlet_bc,
neumann=neumann_bc,
robin=robin_bc,
)
# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
problem.solve()
(2, True)
# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, 0.2))
2022-03-22 21:34:36,870 [754] INFO pulse.iterate: Iterating to target control (f_26)...
2022-03-22 21:34:36,874 [754] INFO pulse.iterate: Current control: f_26 = 0.000,0.000
2022-03-22 21:34:36,874 [754] INFO pulse.iterate: Target: 1.000,0.200
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_DTOL.
*** Warning: PETSc SNES solver diverged in 2 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 2 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 76),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 244),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 377),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 436),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 495),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 554),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 621),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 738),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 797),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 864),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 973),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1074),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1233),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1284),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1359),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1469),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1537)],
[(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 72),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 74)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 240),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 242)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 373),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 375)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 432),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 434)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 491),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 493)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 550),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 552)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 617),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 619)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 734),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 736)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 793),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 795)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 860),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 862)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 969),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 971)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1070),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1072)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1229),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1231)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1280),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1282)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1355),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1357)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1465),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1467)),
(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1533),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 0), FiniteElement('Real', tetrahedron, 0)), 1535))])
# Get the solution
u, p = problem.state.split(deepcopy=True)
volume = geometry.cavity_volume(u=u)
print(f"Cavity volume = {volume}")
Cavity volume = 1.1409709245419588
# Move mesh accoring to displacement
u_int = interpolate(u, dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1))
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, show=False))
fig.show()